Permeability Estimates of Self-AfRne Fracture Faults Based on Generalization of the 

Bottle Neck Concept 



O 
< 



O 



> 



O 

o 



X 



Laurent TalorQ and Harold Aurador(l| 

Univ. Pierre et Marie Curie-Paris6, Univ. Paris-Sud, CNRS, 
Lab. FAST, Bat. 502, Campus Univ., Orsay, F-91405, France. 

Alex Hanser0 

Department of Physics, Norwegian University of Science and Technology, N-l^Ol Trondheim, Norway 

(Dated: April 8, 2010) 

We propose a method for calculating the effective permeability of two-dimensional self-afRne 
permeability fields based on generalizing the one-dimensional concept of a bottleneck. We test the 
method on fracture faults where the local permeability field is given by the cube of the aperture 
field. The method remains accurate even when there is substantial mechanical overlap between the 
two fracture surfaces. The computational efficiency of the method is comparable to calculating a 
simple average and is more than two orders of magnitude faster than solving the Reynolds equations 
using a finite-difference scheme. 

PACS numbers: 



In many low permeability geological formations, flow 
occurs primarily through fracture networks ,19J. In or- 
der to model such systems and to predict their behavior, 
there is a need for reliable modeling of the hydromechan- 
ical behavior of fracture. We consider in this note the 
situation where the shear displacement between the frac- 
ture walls strongly affects its permeability. Because of 
its relevance, this situation has been considered in many 
recent hydromechanical studies 0, i, [13 [H, US [13, [13 . 
Laboratory tests report that the shearing process results 
in a significant channelization of the flow and an enhance- 
ment of the permeability in the direction normal to the 
shear. This behavior is found to be related to the long- 
range spatial organization of the void space, and efforts 
have been undertaken to modelize such system in order to 
provide upscaled value for the fracture permeability. Re- 
cently Mallikanas and Rajaram [iJl determined analyt- 
ically using perturbation analysis of the Reynolds equa- 
tion to the lowest non-trivial order the fracture perme- 
ability. This model, however, does not take into account 
the role of contact areas and will fail if they appear. The 
effect of contacts may, however, be taken into account by 
introducing an empirical parameter [29^ that is strongly 
influenced by the number and the spatial distribution of 
the contacts [T3| . 

We present in this note a computational method for 
calculating the permeability of such fracture faults even 
in the presence of contacts. This new method scales lin- 
early with the number of grid points and is more than 
two orders of magnitude faster than solving the finite- 
differenced Reynolds equations through LU decomposi- 
tion. 
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There is now ample experimental and observational 
evidence that fracture surfaces are self affine, see e.g., 
d [3, H^, [ii, mi, [H, 113. A self-affine fracture may char- 
acterized by a rescaling r of distances in the average frac- 
ture plane and a rescaling of distances in the orthog- 
onal direction leaves the statistical properties of the sur- 
face unchanged. Here C is the Hurst exponent. We use 
in the following C, = 0.8. i.e. the value often reported 
for rocks fractures [l^, [12, [13] ■ When the two matching 
fracture surfaces are displaced by a distance A along the 
average fracture plane, the ensuing aperture field will be 
self-affine up to the length scales of the order of A. On 
larger scales, the aperture field settles to a constant value 
proportional to the average fracture opening j21|. This 
gives rise to a aperture field h = h(x,y). 

As the two fracture surfaces approach each other, they 
will eventually come into contact and hence overlap. 
Overlap also occurs if the gap between the two fracture 
surfaces remains fixed while the lateral displacement A 
increases. At such places of contact, we set the aperture 
h(x,y) to zero. The contact areas are shown as white in 
Fig- HI The lateral displacement results in a strong struc- 
tural anisotropy: The contact zones are more elongated 
in the direction normal to the displacement (x direction 
with reference to Fig. [1]) than in the other direction {y 
direction with reference to Fig. 

The idea behind the method we introduce in this note 
is based on the generalization of the concept of the bot- 
tle neck to higher dimensions. Some forty years ago, [l| 
presented a different generalization of the same concept. 
As we shall see, only in a limiting case do the two gener- 
alizations approach each other. 

Before delving into the two-dimensional generalization 
— and hence the method we present — we discuss the 
one-dimensional case. Hence we have a channel with self- 
affine walls that have been translated relative to each 
other along the direction of the channel by a distance A. 
The channel aperture is given hy h = h{y), where the 
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Figure 1: Aperture field obtained by shifting laterally by A = 
lOe^ two matching self-affine surfaces with Hurst exponent 
C — 0.8. The size is 512 x 512 and the mean aperture is 
(h) — 7. Darker shades mean smaller aperture whereas lighter 
shades means larger apertures. White zones are contact areas. 
Top (resp. bottom) figure: the fiow is normal to (resp. along) 
the lateral displacement 5. The fiow lines are shown as grey 
paths. The worst paths normal to the average fiow directions 
in the two cases are shown as thick grey lines. In the top 
figure, we have (h^)/ d£ ■ eyh{£)^ — 3.32 and in the bottom 
figure {h'^)/ d£- e^hilf = 8.54. 



y axis is oriented along the channel. Assuming that the 
Reynolds equations govern the flow in the channel, the 
local permeability is proportional to h{y)^ . The perme- 
ability of the entire channel is then given by the harmonic 
mean of the local permeability, i.e., oc {h{y)~^)~^ [1^. If 
we now assume that the two channel walls are brought 
close together (so that {h{y)) decreases), the permeabil- 



ity is increasingly controlled by the region of minimum 
aperture imnyh{y)^ [2^, which may be then viewed 
as a bottle neck. 

How wide. A, is the bottle neck region? This of course 
depends on the geometry of the two channel walls in this 
region. For the time being, we leave A as a parameter. 
We now divide the entire channel along the y axis into 
two regions: The bottle neck region which has a width 
A and the rest which has a width L — A, where L is 
the length of the entire channel. The bottle neck region 
have a permeability essentially given by mhiyh{y)^/A 
and the rest of the channel will have a permeability that 
is essentially {h^{y))/{L — A). The total permeability of 
the channel Ky may then be approximated by Dy given 
by 



Dy minyhiyr (/^(y)^) ' " 

Clearly, A will evolve as the average channel width {h{y)) 
decreases and keeping it constant will constitute an ap- 
proximation. How good is such an approximation? A 
natural choice for a fixed A may e.g. be the discretiza- 
tion length scale (i.e., the lattice constant). As {h{y)) 
decreases, the more dominant the bottle neck region will 
be and the more sensitive Dy will be to the discretiza- 
tion at this point. Approximations are unavoidable as 
the average channel width decreases. A "natural" choice 
as the discretization length itself breaks down when the 
discretization itself breaks down. 

We now turn to generalizing this discussion to two- 
dimensional aperture fields h = h{x, y). The main differ- 
ence between the one-dimensional channel and the two- 
dimensional fracture is that flow can in the latter case 
easily bypass regions of small aperture. They do not play 
the crucial role here as they did in the one-dimensional 
channel. 

We therefore generalize the concept of the bottle neck 
for two-dimensional fractures. As a first step to this gen- 
eralization, we consider paths going from one side of the 
fracture to the opposite side cutting across the average 
fiow direction. As a result of mass conservation, the flow 
has to pass through all such paths. For each transverse 
path C with respect to the fiow direction (here, the y 
direction), we may calculate the average aperture cubed 
along it, 

Lc{h')c^ j^M-eMA\ (2) 

where Lq is the length of the path and Cx is the unit 
vector in the x direction. This average now replaces for 
the two dimensional system, the local permeability (y) 
for the channel in one dimension. 

In the one-dimensional channel we then went on to 
identifying the smallest local permeability miuj, /i'^(y). 
This was the bottle neck. In two dimensions, we now 
search for the path with the smallest average aperture 
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cubed, henceforth refered to as the worst path, 



min h{y)" 
y 
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(3) 



Fig. [T] shows for one of the realizations the two worst 
paths obtained for flow directions along and normal to 
the lateral displacement. 

Using the same reasoning as in one dimension, we 
may now generalize Eq. ([1]) by replacing min^ h{y)^ by 
exh{£)^, hence 
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where W is the width of the fracture in the x direction. 

When the average flow is in the x direction rather than 
the y direction, there is of course an equivalent expression 
for Dx- These two expressions, for Dy and D^, form the 
core of our method for approximating the permeability 
of fracture faults. 

We now discuss briefly the relation between the worst 
path method and the Ambegaokar-Halperin-Langer esti- 
mate The AHL estimate is based on the idea that 
when the local permeability is very widely distributed, 
the upscaled permeability is controlled by the smallest 
local permeability along the path connecting the inlet 
to the outlet that has the highest average permeability 
along it. If in Eq. (0), we assume that h^{x, y) is so widely 
distributed that the integral is dominated by the largest 
value of h^{x,y) along the path, the integral becomes 



{h^)c = maxh{i)^ 
tec 



(5) 



If we now combine this expression with Eq. ([3]) to esti- 
mate the permeability of the bottle neck region, we find 
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This expression is essentially the Ambegaokar-Halperin- 
Langer expression for the permeability, except that we 
in this limit end up with the maximum permeability 
along the path with the minimum permeability along it, 
whereas in the analysis of [1], "min" and "max" have 
been substituted. In two-dimensional systems, this is 
equivalent. Hence, only in the limit of extremely broad 
aperture distributions, is our formulation equivalent to 
that of 1]. 

As in one dimension, the width of the bottle neck re- 
gion, A, is a parameter depending on the local topog- 
raphy near the worst path. It needs to the determined 
independently. One way to estimate it is to equate Dy 
(resp. Dx), gotten from Eq. with the permeability 
gotten from another method when the fracture opening 
is large: the detail of the procedure is described farther 
in the text. As in one dimension, we expect A to change 
as the average fracture aperture, {h{x, y)) is lowered. As- 
suming that it is a constant — as we will do — constitutes 
an approximation. 



Given the aperture fields, we compared their perme- 
abilities found using Eq. ^ with the results of two other 
techniques. The first one, proposed by 0, is based on a 
stochastic continuum theory applied to a first order per- 
turbation expansion of the Darcy's law. We compute nu- 
merically the Fourier transform of the permeability field 
perturbation K(kx,ky) = FT{h^{x,y) - {h{x,yf)). The 
two component of the effective permeability are then cal- 
culated from the integrals 
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Since this is only a second order expansion, these results 
are expected to be valid only for small permeability fluc- 
tuations, i.e., when the fracture opening {h) is large com- 
pared to the height fluctuations in the fracture [IJ] . The 
second method consists in solving the flow field inside the 
permeability field by using a lattice Boltzmann method. 
In the this scheme, we introduce a body force to produce 
a Darcy-Brinkman equation as described in [13, . We 
decrease the Brinkman term so that it has no apprecia- 
ble effect on the permeability. When the two surfaces are 
in contact, the lattice site is set to be solid by using the 
"bounce-back" reflection method for the density distri- 
bution. A pressure-imposed boundary condition is used 
at the inlet and outlet as described by [soj . 

Practically, we identify the worst path and the cor- 
responding integral, Eq. (jS]) by using a transfer matrix 
algorithm 4] . If the average flow direction is in the y di- 
rection, the path we construct runs between the sides of 
the sample parallel to the x axis. We discretize the aper- 
ture field h{x,y) — ?> h{i,j), onto a square lattice where 
i runs from 1 to M — W/a and j from 1 to = L/a, 
and a is the lattice constant. We introduce a second field 
p{i,j) which initially is set to zero everywhere. We then 
update layer by layer in the i direction 

p{-i + 1,.?) = 

min[p(i,j - l),p{i,j),p{i,j + 1)] + h{t + , (9) 

until i — M — 1. The integral Eq. ^ is then given by 



P{M,jM) 



min p{M, j) 
j 



(10) 



where we designate by jm the j value where the mininum 
p was identified. In order to reconstruct the worst path, 
we start at the position {M,jM)- We then move on the 
next layer, and identify min[p(M — 1, — — 
l,jM),p{M — l,jM + !)]■ The j value that corresponds 
to the minimum at level i = M — 1 is designated Jm-i- 
We then repeat this algorithm until we have identified ji . 
The sequence ji where i = 1, • • • ,M gives the coordinates 
of the worst path. 
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Figure 2: Normalized permeabilities of a rough fracture for 
flow along the lateral displacement as function of the fracture 
opening {h) for two lateral displacements (Thick solid curves 
A = lOe^, solid curves A = 20e^). Circles: permeability, 
Dy, calculated using Eq. Q. Squares and triangles are for 
the permeabilities obtained by the Lattice Boltzmann [Ky) 
algorithm and by the second-order perturbation theory {Fy). 
The system size is 512 x 512 and A has been set equal to 33 
for A = 10 and 38 for A = 20 by matching Dy and Fy for the 
maximum fracture opening. 
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Figure 3: Normalized permeabilities of a rough fracture for 
flow normal to A as function of the fracture opening (/i) for 
two lateral displacements (Thick solid lines A = lOe^ - Solid 
lines A = 20e^). Despite the flow direction, the conditions are 
similar to the ones of Fig. [2]and circles, squares and triangles 
referred to Dx, and Fx- 

In Eq. ^ we are assuming that the paths only connect 
nearest-neighbor and next-nearest-neighbor nodes on the 
lattice, i.e., (i, j ± fc) with (i + 1, j) where fc = 0, 1. This 
may be generahzed to fc = 0, • • • , m. In our numerical 
calculations presented in Figs.[5]andnwe have used m = 2. 
However, we see no appreciable difference between this 
value and m = 1. 

The algorithm described in Eq. ^ assumes that the 
paths do not turn back, i.e., jc = jc{i)- In very strongly 
disordered fractures, such turns may play a role. This 
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Figure 4: The relative errors \Fx - Kx\/{h^), \Fy - Ky\/{h^), 
\Dx — Kx\/{h^), and \Dy — Ky\/{h^) as a function of the av- 
erage fracture aperture (h). The data has been averaged over 
10 samples, each of size 512 x 512, Hurst exponent H = 0.8 
and lateral displacement u — 10. 



is not the case for the fractures studied here. However, 
when turns do appear, different and more involved al- 
gorithms must be used [ill, ■ Whereas the algorithm 
described in Eq. ^ scales as the number of nodes Af x N 
in the discretized height field, the algorithms capable of 
handling overhangs scales as AP x iV^. 

We first study the situation where flow is parallel to 
the lateral displacement, i.e. orthogonal to the channel- 
ization. Such flow situation is illustrated in top figure in 
Fig. [TJ Fig. [2] shows the variation of the permeability of 
the fracture estimated by the lattice-Boltzmann method 
(squares) and by the second order expansion (triangles) 
as functions of the mean fracture opening and for two 
lateral displacements. For each of the aperture fields, 
p{M,jM) as well as (h^) were measured. The estima- 
tion of Dy still requires an estimation of A, see Eq. (U). 
For the two lateral displacements and for large fracture 
openings, the second order estimate of the permeability, 
Fy, Eq. ^ fits the lattice-Boltzmann Ky well. In this 
region, we equate Fy and Dy, hence determining A. We 
then go on to using the same A for all subsequent frac- 
ture openings. Herein lies the major approximation in 
our method. As soon as contact areas appear (here a no- 
ticeable difference occurs when contacts cover about 10% 
of the total fracture area) , the perturbative estimate Fy 
fails to describe the continuous drop of the permeability 
whereas Dy remains very close to Ky. 

For fiow normal to the lateral displacement and, as il- 
lustrated in the bottom figure of Fig. [U we find strong 
channelization and it is markly different from the one 
observed when flow is parallel to A. Fig. □ shows the per- 
meabilities found by the three methods: a drop off of the 
permeability with the fracture closure is observed but 
is less marked than for flow along the shift (See Fig. [5] 
for comparison). As previously mentioned, as soon as 
contacts between the two surfaces occur the pertuba- 
tive method fails to describe the marked permeability 
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decrease observed with the lattice boltzmann method. 
Yet the worst path method stih accurately captures the 
permeability reduction estimate in the direction parallel 
to the channelization even for large lateral displacement 
of the fracture walls. 

We show in Fig.n the relative errors between the worst 
path method and the lattice Boltzmann method, and the 
perturbative approximation and the lattice Boltzmann 
method for different average fracture openings. The data 
have been averaged over ten samples. As we see, the 
worst path method performs very well for all values of the 
average fracture opening and for the two flow directions. 

To conclude, we have introduced a new technique to 
estimate the permeability of self-affine fracture faults. 
Compared to other approximative methods, it performs 
very well by being able to reproduce the permeability 
closely even when the fracture opening tends to zero. 
"Exact" methods such as the lattice Boltzmann method 
gives more precise results. However, the computation 
time is reduced by several orders of magnitude compared 
to alternative methods. To our knowledge, solving the 
finitc-diffcrcnccd Reynolds equations through LU decom- 
position is the fastest "exact method" . For the samples 
studied in this note, the worst path method used 0.01 
seconds per sample and per average fracture opening, 
whereas the LU decomposition used from 3 to 8 seconds. 
Both methods scale linearly with the number of nodes. 

A length scale A is introduced in order to fit the perme- 
ability measured. This length characterizes the extension 



in the flow direction of the region dominated by the worst 
path. We assume that A remains constant as the aver- 
age fracture opening is changed. This is one of the major 
approximation build into the method — but it allows us 
to determine A by comparison with other approximate 
methods such as the perturative scheme for large enough 
average fracture openings for them to be accurate. 

Future work will be devoted to the study the relation- 
ship of between A and the statistics of the aperture fields. 

The worst path method is accurate even if the aper- 
ture field shows structural anisotropy. Such situation is 
achieved by laterally displacing the fracture walls lead- 
ing, as observed on natural fractures, to an anisotropic 
permeability field. The proposed method can be further 
extended to other transport properties such as diffusion 
or electrical conductivity. Different statistical fields such 
as log normal permeability fields which also give rise to 
heterogenous flow structures will also be of interest as 
well as three-dimensional permeability fields. 
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